# require pacman
if (!require("pacman")) install.packages("pacman")
## Loading required package: pacman
pacman::p_load(
ggplot2,
mvtnorm,
GGally,
corrplot,
readxl,
tidyverse,
gridExtra,
grid,
plotly,
ggcorrplot,
FactoMineR,
factoextra,
rgl,
scatterplot3d,
inspectdf,
mvnormtest,
pastecs,
viridis,
psych
)
data = read_csv("Data/2017_Yellow_Taxi_Trip_Data.csv")
## New names:
## Rows: 22699 Columns: 18
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (3): tpep_pickup_datetime, tpep_dropoff_datetime, store_and_fwd_flag dbl (15):
## ...1, VendorID, passenger_count, trip_distance, RatecodeID, PULoca...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
data = data[, -c(1)]
head(data)
## # A tibble: 6 × 17
## VendorID tpep_pickup_datetime tpep_dropoff_datetime passenger_count
## <dbl> <chr> <chr> <dbl>
## 1 2 3/25/2017 8:55 3/25/2017 9:09 6
## 2 1 4/11/2017 14:53 4/11/2017 15:19 1
## 3 1 12/15/2017 7:26 12/15/2017 7:34 1
## 4 2 5/7/2017 13:17 5/7/2017 13:48 1
## 5 2 4/15/2017 23:32 4/15/2017 23:49 1
## 6 2 3/25/2017 20:34 3/25/2017 20:42 6
## # ℹ 13 more variables: trip_distance <dbl>, RatecodeID <dbl>,
## # store_and_fwd_flag <chr>, PULocationID <dbl>, DOLocationID <dbl>,
## # payment_type <dbl>, fare_amount <dbl>, extra <dbl>, mta_tax <dbl>,
## # tip_amount <dbl>, tolls_amount <dbl>, improvement_surcharge <dbl>,
## # total_amount <dbl>
str(data)
## tibble [22,699 × 17] (S3: tbl_df/tbl/data.frame)
## $ VendorID : num [1:22699] 2 1 1 2 2 2 2 2 2 1 ...
## $ tpep_pickup_datetime : chr [1:22699] "3/25/2017 8:55" "4/11/2017 14:53" "12/15/2017 7:26" "5/7/2017 13:17" ...
## $ tpep_dropoff_datetime: chr [1:22699] "3/25/2017 9:09" "4/11/2017 15:19" "12/15/2017 7:34" "5/7/2017 13:48" ...
## $ passenger_count : num [1:22699] 6 1 1 1 1 6 1 1 1 1 ...
## $ trip_distance : num [1:22699] 3.34 1.8 1 3.7 4.37 ...
## $ RatecodeID : num [1:22699] 1 1 1 1 1 1 1 1 1 1 ...
## $ store_and_fwd_flag : chr [1:22699] "N" "N" "N" "N" ...
## $ PULocationID : num [1:22699] 100 186 262 188 4 161 79 237 234 239 ...
## $ DOLocationID : num [1:22699] 231 43 236 97 112 236 241 114 249 237 ...
## $ payment_type : num [1:22699] 1 1 1 1 2 1 1 1 2 1 ...
## $ fare_amount : num [1:22699] 13 16 6.5 20.5 16.5 9 47.5 16 9 13 ...
## $ extra : num [1:22699] 0 0 0 0 0.5 0.5 1 1 0 0 ...
## $ mta_tax : num [1:22699] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
## $ tip_amount : num [1:22699] 2.76 4 1.45 6.39 0 2.06 9.86 1.78 0 2.75 ...
## $ tolls_amount : num [1:22699] 0 0 0 0 0 0 0 0 0 0 ...
## $ improvement_surcharge: num [1:22699] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 ...
## $ total_amount : num [1:22699] 16.56 20.8 8.75 27.69 17.8 ...
duplicated_rows = data[duplicated(data),]
duplicated_rows
## # A tibble: 0 × 17
## # ℹ 17 variables: VendorID <dbl>, tpep_pickup_datetime <chr>,
## # tpep_dropoff_datetime <chr>, passenger_count <dbl>, trip_distance <dbl>,
## # RatecodeID <dbl>, store_and_fwd_flag <chr>, PULocationID <dbl>,
## # DOLocationID <dbl>, payment_type <dbl>, fare_amount <dbl>, extra <dbl>,
## # mta_tax <dbl>, tip_amount <dbl>, tolls_amount <dbl>,
## # improvement_surcharge <dbl>, total_amount <dbl>
colSums(is.na(data))
## VendorID tpep_pickup_datetime tpep_dropoff_datetime
## 0 0 0
## passenger_count trip_distance RatecodeID
## 0 0 0
## store_and_fwd_flag PULocationID DOLocationID
## 0 0 0
## payment_type fare_amount extra
## 0 0 0
## mta_tax tip_amount tolls_amount
## 0 0 0
## improvement_surcharge total_amount
## 0 0
data_type_pickup = class(data$tpep_pickup_datetime)
data_type_dropoff = class(data$tpep_dropoff_datetime)
print(paste("Data type of tpep_pickup_datetime: ", data_type_pickup))
## [1] "Data type of tpep_pickup_datetime: character"
print(paste("Data type of tpep_dropoff_datetime: ", data_type_dropoff))
## [1] "Data type of tpep_dropoff_datetime: character"
data$tpep_pickup_datetime <- as.POSIXct(data$tpep_pickup_datetime, format = "%m/%d/%Y %H:%M")
data$tpep_dropoff_datetime <- as.POSIXct(data$tpep_dropoff_datetime, format = "%m/%d/%Y %H:%M")
data_type_pickup = class(data$tpep_pickup_datetime)
data_type_dropoff = class(data$tpep_dropoff_datetime)
print(paste("Data type of tpep_pickup_datetime: ", data_type_pickup))
## [1] "Data type of tpep_pickup_datetime: POSIXct"
## [2] "Data type of tpep_pickup_datetime: POSIXt"
print(paste("Data type of tpep_dropoff_datetime: ", data_type_dropoff))
## [1] "Data type of tpep_dropoff_datetime: POSIXct"
## [2] "Data type of tpep_dropoff_datetime: POSIXt"
Thông tin về tập dữ liệu
- ‘VendorID’ - Mã cho biết nhà cung cấp TPEP đã cung cấp hồ sơ.
(1=‘Creative Mobile Technologies, LLC’, 2=‘VeriFone Inc’)
- ‘tpep_pickup_datetime’ - Ngày và giờ khi đồng hồ được bật.
- ‘tpep_dropoff_datetime’ - Ngày và giờ khi đồng hồ được ngắt.
- ‘passenger_count’ - Số lượng hành khách trên xe.
- ‘trip_distance’ - Khoảng cách chuyến đi đã trôi qua tính bằng dặm do
đồng hồ taxi báo cáo.
- ‘RatecodeID’ - Mã giá cước cuối cùng có hiệu lực vào cuối chuyến đi.
(1=‘Giá chuẩn’, 2=‘JFK’, 3=‘Newark’, 4=‘Nassau of Westchester’, 5=‘Giá
đã thương lượng’, 6=‘Đi theo nhóm’, 99=‘NULL/không xác định’)
- ‘store_and_fwd_flag’ - Cờ này cho biết liệu hồ sơ chuyến đi có được
lưu trong bộ nhớ xe trước khi gửi đến nhà cung cấp hay không, còn gọi là
“lưu trữ và chuyển tiếp”, vì xe không có kết nối với máy chủ. (Y=‘lưu
trữ và chuyển tiếp chuyến đi’, N=‘không phải là chuyến đi lưu trữ và
chuyển tiếp’)
- ‘PULocationID’ - Khu vực taxi TLC mà đồng hồ tính cước đã được
bật
- ‘DOLocationID’ - Khu vực taxi TLC mà đồng hồ tính cước đã được
tắt
- ‘payment_type’ - Mã số biểu thị cách hành khách thanh toán cho chuyến
đi. (0= ‘Flex Fare trip’, 1= ‘Credit card’, 2= ‘Cash’, 3= ‘No charge’,
4= ‘Dispute’, 5= ‘Unknown’, 6= ‘Voided trip’)
- ‘fare_amount’ - Giá vé theo thời gian và khoảng cách được tính theo
đồng hồ. Để biết thêm thông tin về các cột sau
- ‘extra’ - Các khoản phụ phí và phụ thu khác.
- ‘mta_tax’ - Thuế được tự động kích hoạt dựa trên mức giá theo đồng hồ
đang sử dụng.
- ‘tip_amount’ - Số tiền boa – Trường này được tự động điền cho tiền boa
bằng thẻ tín dụng. Tiền boa bằng tiền mặt không được bao gồm.
- ‘tolls_amount’ - Tổng số tiền của tất cả các khoản phí cầu đường đã
thanh toán trong chuyến đi.
- ‘improvement_surcharge’ - Phụ phí cải thiện được đánh giá cho các
chuyến đi tại thời điểm hạ cờ. Phụ phí cải thiện bắt đầu được áp dụng
vào năm 2015.
- ‘total_amount’ - Tổng số tiền tính cho hành khách. Không bao gồm tiền
boa bằng tiền mặt.
Ghi chú về Bộ dữ liệu: Bộ dữ liệu chứa 22699 hàng và 17 cột, không có giá trị bị thiếu hoặc hàng trùng lặp. Tất cả các cột số đều là số dương.
eda = data
eda = eda[, -c(8,9)]
# Convert vendorID 1='Creative Mobile Technologies, LLC', 2='VeriFone Inc'
eda$VendorID <- ifelse(eda$VendorID == 1, 'Creative Mobile Technologies, LLC', 'VeriFone Inc')
# Convert RatecodeID 1='Standard rate', 2='JFK', 3='Newark', 4='Nassau or Westchester', 5='Negotiated fare', 6='Group ride', 99='NULL/Unknown'
eda$RatecodeID <- ifelse(eda$RatecodeID == 1, 'Standard rate', ifelse(eda$RatecodeID == 2, 'JFK', ifelse(eda$RatecodeID == 3, 'Newark', ifelse(eda$RatecodeID == 4, 'Nassau or Westchester', ifelse(eda$RatecodeID == 5, 'Negotiated fare', ifelse(eda$RatecodeID == 6, 'Group ride', 'NULL/Unknown'))))))
# Convert store_and_fwd_flag Y='store and forward trip', N='not a store and forward trip'
eda$store_and_fwd_flag <- ifelse(eda$store_and_fwd_flag == 'Y', 'store and forward trip', 'not a store and forward trip')
# Convert payment_type 0='Flex Fare trip', 1='Credit card', 2='Cash', 3='No charge', 4='Dispute', 5='Unknown', 6='Voided trip'
eda$payment_type <- ifelse(eda$payment_type == 0, 'Flex Fare trip', ifelse(eda$payment_type == 1, 'Credit card', ifelse(eda$payment_type == 2, 'Cash', ifelse(eda$payment_type == 3, 'No charge', ifelse(eda$payment_type == 4, 'Dispute', ifelse(eda$payment_type == 5, 'Unknown', 'Voided trip'))))))
# Convert to factor
eda$VendorID <- as.factor(eda$VendorID)
eda$RatecodeID <- as.factor(eda$RatecodeID)
eda$store_and_fwd_flag <- as.factor(eda$store_and_fwd_flag)
eda$payment_type <- as.factor(eda$payment_type)
str(eda)
## tibble [22,699 × 15] (S3: tbl_df/tbl/data.frame)
## $ VendorID : Factor w/ 2 levels "Creative Mobile Technologies, LLC",..: 2 1 1 2 2 2 2 2 2 1 ...
## $ tpep_pickup_datetime : POSIXct[1:22699], format: "2017-03-25 08:55:00" "2017-04-11 14:53:00" ...
## $ tpep_dropoff_datetime: POSIXct[1:22699], format: "2017-03-25 09:09:00" "2017-04-11 15:19:00" ...
## $ passenger_count : num [1:22699] 6 1 1 1 1 6 1 1 1 1 ...
## $ trip_distance : num [1:22699] 3.34 1.8 1 3.7 4.37 ...
## $ RatecodeID : Factor w/ 6 levels "JFK","Nassau or Westchester",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ store_and_fwd_flag : Factor w/ 2 levels "not a store and forward trip",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ payment_type : Factor w/ 4 levels "Cash","Credit card",..: 2 2 2 2 1 2 2 2 1 2 ...
## $ fare_amount : num [1:22699] 13 16 6.5 20.5 16.5 9 47.5 16 9 13 ...
## $ extra : num [1:22699] 0 0 0 0 0.5 0.5 1 1 0 0 ...
## $ mta_tax : num [1:22699] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
## $ tip_amount : num [1:22699] 2.76 4 1.45 6.39 0 2.06 9.86 1.78 0 2.75 ...
## $ tolls_amount : num [1:22699] 0 0 0 0 0 0 0 0 0 0 ...
## $ improvement_surcharge: num [1:22699] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 ...
## $ total_amount : num [1:22699] 16.56 20.8 8.75 27.69 17.8 ...
# Extract number columns
num_cols <- eda %>% select_if(is.numeric) %>% names()
num_cols
## [1] "passenger_count" "trip_distance" "fare_amount"
## [4] "extra" "mta_tax" "tip_amount"
## [7] "tolls_amount" "improvement_surcharge" "total_amount"
# stat.desc
stats = stat.desc(eda[,num_cols])
formatted_stats <- format(stats, scientific = FALSE)
formatted_stats
## passenger_count trip_distance fare_amount extra
## nbr.val 22699.000000000 22699.00000000 22699.00000000 22699.000000000
## nbr.null 33.000000000 148.00000000 6.00000000 11921.000000000
## nbr.na 0.000000000 0.00000000 0.00000000 0.000000000
## min 0.000000000 0.00000000 -120.00000000 -1.000000000
## max 6.000000000 33.96000000 999.99000000 4.500000000
## range 6.000000000 33.96000000 1119.99000000 5.500000000
## sum 37279.000000000 66129.29000000 295691.46000000 7565.000000000
## median 1.000000000 1.61000000 9.50000000 0.000000000
## mean 1.642319045 2.91331292 13.02662937 0.333274594
## SE.mean 0.008530566 0.02424748 0.08790406 0.003073748
## CI.mean.0.95 0.016720495 0.04752673 0.17229798 0.006024756
## var 1.651819029 13.34565969 175.39798725 0.214458441
## std.dev 1.285231119 3.65317118 13.24379052 0.463096579
## coef.var 0.782570916 1.25395770 1.01667056 1.389534599
## mta_tax tip_amount tolls_amount
## nbr.val 22699.0000000000 22699.00000000 22699.00000000
## nbr.null 90.0000000000 8057.00000000 21525.00000000
## nbr.na 0.0000000000 0.00000000 0.00000000
## min -0.5000000000 0.00000000 0.00000000
## max 0.5000000000 200.00000000 19.10000000
## range 1.0000000000 200.00000000 19.10000000
## sum 11291.5000000000 41670.40000000 7094.38000000
## median 0.5000000000 1.35000000 0.00000000
## mean 0.4974448214 1.83578131 0.31254152
## SE.mean 0.0002619441 0.01858882 0.00928710
## CI.mean.0.95 0.0005134284 0.03643536 0.01820335
## var 0.0015574852 7.84350752 1.95779403
## std.dev 0.0394649873 2.80062627 1.39921193
## coef.var 0.0793354069 1.52557729 4.47688334
## improvement_surcharge total_amount
## nbr.val 22699.0000000000 22699.0000000
## nbr.null 6.0000000000 4.0000000
## nbr.na 0.0000000000 0.0000000
## min -0.3000000000 -120.3000000
## max 0.3000000000 1200.2900000
## range 0.6000000000 1320.5900000
## sum 6799.5000000000 370232.0900000
## median 0.3000000000 11.8000000
## mean 0.2995506410 16.3105022
## SE.mean 0.0001040259 0.1068439
## CI.mean.0.95 0.0002038979 0.2094213
## var 0.0002456347 259.1229160
## std.dev 0.0156727376 16.0972953
## coef.var 0.0523208283 0.9869282
Có một số thông tin nổi bật từ bảng thống kê tóm tắt này. Rõ ràng có
một số giá trị ngoại lệ trong một số biến, như tip_amount
($200) và total_amount ($1.200). Ngoài ra, một số biến, như
mta_tax, dường như gần như không đổi trong toàn bộ dữ liệu,
ta không mong đợi chúng mang lại khả năng dự đoán cao.
plot_all_densities <- function(data, include_cols = num_cols , fill_color = "skyblue", alpha_value = 0.3, base_size = 15) {
plot_list <- list()
for (var in names(data[, num_cols])) {
p = ggplot(data, aes_string(x = var)) +
geom_histogram(fill = "skyblue", color = "blue", alpha = alpha_value, stat='count') +
theme_minimal(base_size = base_size) +
labs(title = var,
x = "",
y = "Frequency",
) +
geom_rug(sides = "b")
plot_list[[var]] <- p
}
title_grob <- textGrob("Distribution Of Numerical Columns", gp = gpar(fontsize = 20, fontface = "bold"))
grid.arrange(title_grob, grobs=plot_list, ncol = 3, nrow = 3)
}
plot_all_densities(eda)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning in geom_histogram(fill = "skyblue", color = "blue", alpha = alpha_value, : Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
## Ignoring unknown parameters: `binwidth`, `bins`, and `pad`
# countplot for factor columns
plot_all_countplots <- function(data, include_cols = c("VendorID", "RatecodeID", "store_and_fwd_flag", "payment_type"), fill_color = "skyblue", alpha_value = 0.8, base_size = 15) {
plot_list <- list()
for (var in include_cols) {
p = ggplot(data, aes_string(x = var)) +
geom_bar(fill = "aquamarine4", alpha = alpha_value) +
theme_minimal(base_size = base_size) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(
plot.title = element_text(size = 16),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16),
legend.title = element_text(size = 16),
legend.text = element_text(size = 16)
) +
labs(title = var,
x = "",
y = "Frequency",
) +
geom_text(stat='count', aes(label=..count..), vjust=-0.5)
plot_list[[var]] <- p
}
title_grob <- textGrob("Countplot For Factor Columns", gp = gpar(fontsize = 20, fontface = "bold"))
grid.arrange(title_grob, grobs=plot_list, ncol = 2, nrow = 2)
}
plot_all_countplots(eda)
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
duration, thể hiện tổng số phút
của mỗi chuyến đi taxi.eda$duration <- as.numeric(difftime(eda$tpep_dropoff_datetime, eda$tpep_pickup_datetime, units = "mins"))
head(eda$duration)
## [1] 14 26 8 31 17 8
trip_distance,
fare_amount, durationnum_cols <- eda %>% select_if(is.numeric) %>% names()
num_cols
## [1] "passenger_count" "trip_distance" "fare_amount"
## [4] "extra" "mta_tax" "tip_amount"
## [7] "tolls_amount" "improvement_surcharge" "total_amount"
## [10] "duration"
boxplot(eda[, num_cols], rotate = 45, col = "skyblue", main = "Boxplot of numerical columns")
trip_distance outliershead(sort(unique(eda$trip_distance)), 10)
## [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
Khoảng cách được ghi lại với độ chính xác cao. Tuy nhiên, có thể các chuyến đi có khoảng cách bằng không nếu hành khách gọi taxi rồi đổi ý. Bên cạnh đó, có đủ giá trị bằng không trong dữ liệu để gây ra vấn đề không?
sum(eda$trip_distance == 0)
## [1] 148
148 trong số ~23.000 chuyến đi là không đáng kể. Ta có thể quy cho nó
giá trị 0,01, nhưng nó không có nhiều tác động đến mô hình. Do đó, cột
trip_distance sẽ không bị ảnh hưởng đối với các giá trị
ngoại lệ.
fare_amountsummary(eda$fare_amount)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -120.00 6.50 9.50 13.03 14.50 999.99
Phạm vi giá trị trong cột fare_amount rất lớn và các giá
trị cực trị cũng không có nhiều ý nghĩa.
* Với các giá trị thấp: Giá trị âm là có vấn đề. Giá trị bằng 0 có thể
hợp lệ nếu taxi ghi lại chuyến đi đã bị hủy ngay lập tức.
* Với các giá trị cao: Số tiền cước tối đa trong tập dữ liệu này gần
1.000 đô la, điều này có vẻ rất khó xảy ra trong một chuyến taxi. Giá
trị cao cho tính năng này có thể được giới hạn dựa trên trực giác và số
liệu thống kê. Phạm vi tứ phân vị (IQR) là 8 đô la. Công thức chuẩn Q3 +
(1,5 * IQR) cho ra kết quả là 26,50 đô la. Có vẻ không phù hợp với mức
giới hạn cước tối đa. Trong trường hợp này, chúng ta sẽ sử dụng hệ số 6
* IQR, với kết quả ở mức 6*IQR là 62,5 đô la
Với các giá trị nhỏ hơn 0: Ta quy về 0
eda$fare_amount[eda$fare_amount < 0] <- 0
min(eda$fare_amount)
## [1] 0
Tiếp theo gán giá trị lớn nhất là Q3 + (6 * IQR).
# Tiếp theo gán giá trị lớn nhất là `Q3 + (6 * IQR)`.
Q3 <- quantile(eda$fare_amount, 0.75)
IQR <- IQR(eda$fare_amount)
upper_bound <- Q3 + (6 * IQR)
eda$fare_amount[eda$fare_amount > upper_bound] <- upper_bound
summary(eda$fare_amount)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 6.5 9.5 12.9 14.5 62.5
durationsummary(eda$duration)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -17.00 7.00 11.00 17.01 18.00 1440.00
eda$duration[eda$duration < 0] <- 0
min(eda$duration)
## [1] 0
Q3 <- quantile(eda$duration, 0.75)
IQR <- IQR(eda$duration)
upper_bound <- Q3 + (6 * IQR)
eda$duration[eda$duration > upper_bound] <- upper_bound
summary(eda$duration)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 7.00 11.00 14.44 18.00 84.00
total_amountsummary(eda$total_amount)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -120.30 8.75 11.80 16.31 17.80 1200.29
Cột total_amount cũng có các giá trị ngoại lệ ở cả hai
extreme value dưới và trên.
* Giá trị thấp: Giá trị âm không hợp lý. Ta sẽ gán tất cả giá trị âm
bằng 0.
* Giá trị cao: Gán các giá trị cao theo cùng phương pháp đã gán các giá
trị ngoại lệ cho giá vé: Q3 + (6 * IQR).
eda$total_amount[eda$total_amount < 0] <- 0
min(eda$total_amount)
## [1] 0
Q3 <- quantile(eda$total_amount, 0.75)
IQR <- IQR(eda$total_amount)
upper_bound <- Q3 + (6 * IQR)
eda$total_amount[eda$total_amount > upper_bound] <- upper_bound
summary(eda$total_amount)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 8.75 11.80 16.12 17.80 72.10
boxplot(eda[, num_cols], rotate = 45, col = "skyblue", main = "Boxplot of numerical columns")
eda$day_of_month <- as.numeric(format(eda$tpep_pickup_datetime, "%d"))
taxi_usage <- table(eda$day_of_month)
taxi_usage <- as.data.frame(taxi_usage)
colnames(taxi_usage) <- c("Day", "Count")
ggplot(taxi_usage, aes(x = Day, y = Count, group=1)) +
geom_point(color = "coral4") +
geom_line(color = "coral3") +
labs(
title = "Figure 1: How does taxi usage change over a month",
x = "Day of the Month",
y = "Number of Taxi Rides"
) +
theme_minimal()
Chúng tôi quan sát thấy một mô hình sử dụng taxi tương đối ổn định trong suốt mỗi ngày trong tháng. Tuy nhiên, có sự sụt giảm đáng kể về số chuyến đi vào cuối tháng. Có vẻ như họ đã sử dụng hết ngân sách đi lại của mình vào đầu và giữa tháng.
eda$month <- as.numeric(format(eda$tpep_pickup_datetime, "%m"))
taxi_usage_month <- table(eda$month)
taxi_usage_month <- as.data.frame(taxi_usage_month)
colnames(taxi_usage_month) <- c("Month", "Count")
ggplot(taxi_usage_month, aes(x = factor(Month), y = Count, fill = factor(Month))) +
geom_bar(stat = "identity") +
geom_text(aes(label = Count), vjust = -0.5, color = "black") +
scale_fill_viridis_d() +
labs(
title = "Figure 2: Which month has the most taxi rides in the year",
x = "Month",
y = "Number of Taxi Rides"
) +
theme_minimal() +
theme(legend.position = "none")
taxi_usage_vendor <- table(eda$VendorID)
taxi_usage_vendor <- as.data.frame(taxi_usage_vendor)
colnames(taxi_usage_vendor) <- c("Vendor", "Count")
ggplot(taxi_usage_vendor, aes(x = Vendor, y = Count, fill = Vendor)) +
geom_bar(stat = "identity") +
geom_text(aes(label = Count), vjust = -0.5, color = "black") +
scale_fill_viridis_d() +
labs(
title = "Figure 3: Which vendor has the most taxi rides",
x = "Vendor",
y = "Number of Taxi Rides"
) +
theme_minimal() +
theme(legend.position = "none")
Nhận xét: Có sự khác biệt về khoảng cách chuyến đi giữa các nhà cung cấp, nhưng không quá lớn.
total_amountgrid.arrange(
ggplot(eda, aes(x = VendorID, y = trip_distance, fill = VendorID)) +
geom_boxplot() +
labs(
title = "Figure 4: Is there a difference in trip distance between vendors",
x = "Vendor",
y = "Trip Distance"
) +
theme_minimal(),
ggplot(eda, aes(x = RatecodeID, y = total_amount, fill = RatecodeID)) +
geom_boxplot() +
labs(
title = "Figure 5: Influence of RatecodeID on total amount",
x = "RatecodeID",
y = "Total Amount"
) +
theme_minimal(),
ggplot(eda, aes(x = store_and_fwd_flag, y = total_amount, fill = store_and_fwd_flag)) +
geom_boxplot() +
labs(
title = "Figure 6: Influence of store_and_fwd_flag on total amount",
x = "store_and_fwd_flag",
y = "Total Amount"
) +
theme_minimal(),
ggplot(eda, aes(x = payment_type, y = total_amount, fill = payment_type)) +
geom_boxplot() +
labs(
title = "Figure 7: Influence of payment_type on total amount",
x = "payment_type",
y = "Total Amount"
) +
theme_minimal(),
ncol = 2
)
eda$day_of_week <- as.numeric(format(eda$tpep_pickup_datetime, "%u"))
eda$quarter <- as.numeric(format(eda$tpep_pickup_datetime, "%m")) %/% 4 + 1
daily_revenue <- eda %>%
group_by(day_of_week) %>%
summarise(total_amount = sum(fare_amount, na.rm = TRUE)) %>%
arrange(day_of_week)
monthly_revenue <- eda %>%
group_by(month) %>%
summarise(total_amount = sum(fare_amount, na.rm = TRUE)) %>%
arrange(month)
quarterly_revenue <- eda %>%
group_by(quarter) %>%
summarise(total_amount = sum(fare_amount, na.rm = TRUE)) %>%
arrange(quarter)
highest_daily_revenue <- daily_revenue %>%
filter(total_amount == max(total_amount, na.rm = TRUE))
highest_monthly_revenue <- monthly_revenue %>%
filter(total_amount == max(total_amount, na.rm = TRUE))
highest_quarterly_revenue <- quarterly_revenue %>%
filter(total_amount == max(total_amount, na.rm = TRUE))
print(paste("Ngày có doanh thu cao nhất trong tuần là thứ 5",
".Với số tiền là:", highest_daily_revenue$total_amount))
## [1] "Ngày có doanh thu cao nhất trong tuần là thứ 5 .Với số tiền là: 45180.6"
print(paste("Tháng có doanh thu cao nhất trong năm: ", highest_monthly_revenue$month,
"Với số tiền là: ", highest_monthly_revenue$total_amount))
## [1] "Tháng có doanh thu cao nhất trong năm: 5 Với số tiền là: 26823.01"
print(paste("Qúy có doanh thu cao nhất trong năm là quý", highest_quarterly_revenue$quarter,
"Với số tiền là", highest_quarterly_revenue$total_amount))
## [1] "Qúy có doanh thu cao nhất trong năm là quý 2 Với số tiền là 99782.81"
correlation_matrix <- cor(eda[, num_cols])
ggcorrplot(correlation_matrix, lab = TRUE, method = "square", title = "Figure 8: Correlation heatmap of numerical columns")
manova_results <- manova(cbind(trip_distance, fare_amount, duration, total_amount) ~ VendorID, data = eda)
summary(manova_results)
## Df Pillai approx F num Df den Df Pr(>F)
## VendorID 1 0.00037978 2.1555 4 22694 0.07131 .
## Residuals 22697
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Kiểm định ANOVA cho total_amount")
## [1] "Kiểm định ANOVA cho total_amount"
anova_results <- aov(total_amount ~ VendorID, data = eda)
summary(anova_results)
## Df Sum Sq Mean Sq F value Pr(>F)
## VendorID 1 17 16.62 0.1 0.752
## Residuals 22697 3780099 166.55
anova_results <- aov(duration ~ VendorID, data = eda)
print("Kiểm định ANOVA cho duration")
## [1] "Kiểm định ANOVA cho duration"
summary(anova_results)
## Df Sum Sq Mean Sq F value Pr(>F)
## VendorID 1 53 53.27 0.379 0.538
## Residuals 22697 3190306 140.56
anova_results <- aov(fare_amount ~ VendorID, data = eda)
print("Kiểm định ANOVA cho fare_amount")
## [1] "Kiểm định ANOVA cho fare_amount"
summary(anova_results)
## Df Sum Sq Mean Sq F value Pr(>F)
## VendorID 1 3 2.76 0.025 0.875
## Residuals 22697 2522099 111.12
anova_results <- aov(trip_distance ~ VendorID, data = eda)
print("Kiểm định ANOVA cho trip_distance")
## [1] "Kiểm định ANOVA cho trip_distance"
summary(anova_results)
## Df Sum Sq Mean Sq F value Pr(>F)
## VendorID 1 18 18.34 1.375 0.241
## Residuals 22697 302901 13.35
data.pca <- PCA(eda[, num_cols], graph=F)
eig.val <- get_eigenvalue(data.pca)
eig.val
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 4.496876686 44.96876686 44.96877
## Dim.2 1.643460437 16.43460437 61.40337
## Dim.3 1.002979679 10.02979679 71.43317
## Dim.4 0.974999128 9.74999128 81.18316
## Dim.5 0.674774045 6.74774045 87.93090
## Dim.6 0.588084482 5.88084482 93.81174
## Dim.7 0.316920815 3.16920815 96.98095
## Dim.8 0.219676290 2.19676290 99.17772
## Dim.9 0.077175198 0.77175198 99.94947
## Dim.10 0.005053241 0.05053241 100.00000
fviz_eig(data.pca, addlabels = TRUE, ylim = c(0, 100)) +
labs(title = "Scree Plot", x = "Principal Component", y = "Percentage of Variance") +
theme_minimal()
loadings <- as.data.frame(data.pca$var$coord)
loadings$Variable <- rownames(loadings)
loadings
## Dim.1 Dim.2 Dim.3 Dim.4
## passenger_count 0.0157836723 -0.007700298 0.963301835 0.26357361
## trip_distance 0.9315639649 0.061576239 0.015392126 -0.01897864
## fare_amount 0.9664719333 0.039491509 0.016399632 -0.03472738
## extra 0.1393237374 0.156149640 -0.266129046 0.93972532
## mta_tax -0.1762988861 0.893668718 0.018737831 -0.05074638
## tip_amount 0.6606644141 -0.016653396 -0.042364596 -0.05371105
## tolls_amount 0.7269844610 -0.105736408 0.006839232 -0.04751297
## improvement_surcharge 0.0002759055 0.889661305 0.031990056 -0.11326393
## total_amount 0.9854914266 0.043918282 -0.000493820 -0.01014541
## duration 0.8414739681 0.100696906 0.021965234 -0.01513325
## Dim.5 Variable
## passenger_count 0.04359601 passenger_count
## trip_distance -0.15577176 trip_distance
## fare_amount -0.15424734 fare_amount
## extra 0.04069800 extra
## mta_tax -0.10255222 mta_tax
## tip_amount 0.58928165 tip_amount
## tolls_amount 0.30264521 tolls_amount
## improvement_surcharge 0.20360820 improvement_surcharge
## total_amount -0.03550355 total_amount
## duration -0.36204872 duration
p1 <- ggplot(loadings, aes(x = reorder(Variable, Dim.1), y = Dim.1)) +
geom_bar(stat = "identity", fill = "cornsilk4") +
coord_flip() +
labs(title = "Loadings on PC1", x = "Variables", y = "Loadings") +
theme_minimal() +
theme(
plot.title = element_text(size = 14),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14)
)
p2 <- ggplot(loadings, aes(x = reorder(Variable, Dim.2), y = Dim.2)) +
geom_bar(stat = "identity", fill = "coral3") +
coord_flip() +
labs(title = "Loadings on PC2", x = "Variables", y = "Loadings") +
theme_minimal() +
theme(
plot.title = element_text(size = 14),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14)
)
p3 <- ggplot(loadings, aes(x = reorder(Variable, Dim.3), y = Dim.3)) +
geom_bar(stat = "identity", fill = "aquamarine4") +
coord_flip() +
labs(title = "Loadings on PC3", x = "Variables", y = "Loadings") +
theme_minimal() +
theme(
plot.title = element_text(size = 14),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14)
)
grid.arrange(p1, p2, p3, nrow = 2, ncol = 2)
contrib_PC1 <- as.data.frame(data.pca$var$contrib[,1])
colnames(contrib_PC1) <- c("Contribution")
contrib_PC1$Variable <- rownames(contrib_PC1)
contrib_PC2 <- as.data.frame(data.pca$var$contrib[,2])
colnames(contrib_PC2) <- c("Contribution")
contrib_PC2$Variable <- rownames(contrib_PC2)
contrib_PC3 <- as.data.frame(data.pca$var$contrib[,3])
colnames(contrib_PC3) <- c("Contribution")
contrib_PC3$Variable <- rownames(contrib_PC3)
p1 = ggplot(contrib_PC1, aes(x = reorder(Variable, Contribution), y = Contribution)) +
geom_bar(stat = "identity", fill = "skyblue") +
coord_flip() +
labs(title = "Contribution of Variables to Dim 1", x = "", y = "Contribution (%)") +
theme_minimal()
p2 = ggplot(contrib_PC2, aes(x = reorder(Variable, Contribution), y = Contribution)) +
geom_bar(stat = "identity", fill = "salmon") +
coord_flip() +
labs(title = "Contribution of Variables to Dim 2", x = "", y = "Contribution (%)") +
theme_minimal()
p3 = ggplot(contrib_PC3, aes(x = reorder(Variable, Contribution), y = Contribution)) +
geom_bar(stat = "identity", fill = "purple") +
coord_flip() +
labs(title = "Contribution of Variables to Dim 3", x = "", y = "Contribution (%)") +
theme_minimal()
grid.arrange(p1,
p2,
p3,
nrow=1)
p1 <- fviz_pca_var(data.pca,
col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE)
p2 <- fviz_pca_var(data.pca, axes = c(2, 3),
col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE)
grid.arrange(p1, p2, nrow = 1)
fviz_pca_biplot(data.pca, axes = c(1, 2), geom = "point", col.var = "contrib")
fviz_pca_biplot(data.pca, axes = c(2, 3), geom = "point", col.var = "contrib")
# merge biplot pc1 pc2 pc3
biplot_pc1_pc2 <- fviz_pca_biplot(data.pca, axes = c(1, 2), geom = "point", col.var = "contrib")
biplot_pc2_pc3 <- fviz_pca_biplot(data.pca, axes = c(2, 3), geom = "point", col.var = "contrib")
grid.arrange(biplot_pc1_pc2, biplot_pc2_pc3, nrow = 1)
# pca 3d
explained_variance <- data.pca$eig[, 2]
pca_3d <- as.data.frame(data.pca$ind$coord)
pca_3d$VendorID <- data$VendorID
plot_ly(pca_3d, x = ~Dim.1, y = ~Dim.2, z = ~Dim.3, color = ~VendorID, colors = c("#00AFBB", "salmon")) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = paste0('Dim1 (', round(explained_variance[1], 2), '%)')),
yaxis = list(title = paste0('Dim2 (', round(explained_variance[2], 2), '%)')),
zaxis = list(title = paste0('Dim3 (', round(explained_variance[3], 2), '%)'))),
title = "3D PCA Plot of Psychological Data")
eda
## # A tibble: 22,699 × 20
## VendorID tpep_pickup_datetime tpep_dropoff_datetime passenger_count
## <fct> <dttm> <dttm> <dbl>
## 1 VeriFone Inc 2017-03-25 08:55:00 2017-03-25 09:09:00 6
## 2 Creative Mobile T… 2017-04-11 14:53:00 2017-04-11 15:19:00 1
## 3 Creative Mobile T… 2017-12-15 07:26:00 2017-12-15 07:34:00 1
## 4 VeriFone Inc 2017-05-07 13:17:00 2017-05-07 13:48:00 1
## 5 VeriFone Inc 2017-04-15 23:32:00 2017-04-15 23:49:00 1
## 6 VeriFone Inc 2017-03-25 20:34:00 2017-03-25 20:42:00 6
## 7 VeriFone Inc 2017-05-03 19:04:00 2017-05-03 20:03:00 1
## 8 VeriFone Inc 2017-08-15 17:41:00 2017-08-15 18:03:00 1
## 9 VeriFone Inc 2017-02-04 16:17:00 2017-02-04 16:29:00 1
## 10 Creative Mobile T… 2017-11-10 15:20:00 2017-11-10 15:40:00 1
## # ℹ 22,689 more rows
## # ℹ 16 more variables: trip_distance <dbl>, RatecodeID <fct>,
## # store_and_fwd_flag <fct>, payment_type <fct>, fare_amount <dbl>,
## # extra <dbl>, mta_tax <dbl>, tip_amount <dbl>, tolls_amount <dbl>,
## # improvement_surcharge <dbl>, total_amount <dbl>, duration <dbl>,
## # day_of_month <dbl>, month <dbl>, day_of_week <dbl>, quarter <dbl>
Thông thường, trước khi tiến hành phân tích thành phần chính (PCA) và
EFA, ta phải kiểm tra xem dữ liệu ta thu thập được có phù hợp cho phân
tích hay không. Hai phương pháp kiểm định thường dùng để làm điều này là
kiểm định KMO và Bartlett
## Kiểm định KMO
kmo_fa <- KMO(eda[, num_cols])
kmo_fa
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = eda[, num_cols])
## Overall MSA = 0.64
## MSA for each item =
## passenger_count trip_distance fare_amount
## 0.70 0.95 0.61
## extra mta_tax tip_amount
## 0.17 0.33 0.50
## tolls_amount improvement_surcharge total_amount
## 0.60 0.35 0.62
## duration
## 0.96
Một bộ dữ liệu được xem là phù hợp để phân tích EFA thường có giá trị KMO Test tối thiểu là 0.5. Trong trường hợp này, giá trị KMO Test là 0.64, vì vậy ta có thể tiến hành phân tích nhân tố cho bộ dữ liệu.
cortest.bartlett(eda[, num_cols])
## R was not square, finding R from data
## $chisq
## [1] 214695.6
##
## $p.value
## [1] 0
##
## $df
## [1] 45
Bartlett Test có giá trị 19334.49, bậc tự do là 45 và trị số p (p value =0) nhỏ hơn <5%, có nghĩa là quan hệ giữa các biến (items) đủ lớn để sử dụng phân tích EFA.
Trong phân tích EFA, căn cứ để xác định các nhân tố chính được rút ra là sử dụng giá trị của Eigenvalue. Theo tiêu chuẩn của Kaiser thì nhân tố chính được rút ra phải có Eigenvalue > 1. Một tiêu chuẩn khác ít nghiêm ngặt hơn đó là tiêu chuẩn của Jolliffe, theo Jolliffe thì các nhân tố có Eigenvalue > 0.7 có thể được chọn. Trong bộ dữ liệu này, số lượng nhân tố chính được rút ra dựa vào tiêu chuẩn của Kaiser.
fviz_eig(data.pca, addlabels = TRUE, ylim = c(0, 100), n=10) +
labs(title = "Scree Plot", x = "Principal Component", y = "Percentage of Variance") +
theme_minimal()
Xác định các biến cấu thành nhân tố được rút ra, đặt tên nhân tố
Ở trên, sử dụng Eigenvalue theo tiêu chuẩn Kaiser ta đã trích ra được
3 nhân tố chính, để biết các nhân tố này được cấu thành từ những biến
(items) nào ta sử dụng phép xoay Varimax.
print.psych(pc2, cut = 3, sort = TRUE) # Loại bỏ các biến có hệ số tải
(Factor Loading) bé hơn 0.4
pc2 = principal(eda[, num_cols], nfactors = 3, rotate = "varimax")
print.psych(pc2, cut = 0.4, sort = TRUE)
## Principal Components Analysis
## Call: principal(r = eda[, num_cols], nfactors = 3, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
## item RC1 RC2 RC3 h2 u2 com
## total_amount 9 0.98 0.97 0.027 1.0
## fare_amount 3 0.97 0.94 0.064 1.0
## trip_distance 2 0.93 0.87 0.128 1.0
## duration 10 0.84 0.72 0.281 1.0
## tolls_amount 7 0.72 0.54 0.460 1.1
## tip_amount 6 0.66 0.44 0.561 1.0
## mta_tax 5 0.90 0.83 0.170 1.1
## improvement_surcharge 8 0.89 0.79 0.207 1.0
## passenger_count 1 0.96 0.93 0.072 1.0
## extra 4 0.11 0.885 1.8
##
## RC1 RC2 RC3
## SS loadings 4.48 1.64 1.02
## Proportion Var 0.45 0.16 0.10
## Cumulative Var 0.45 0.61 0.71
## Proportion Explained 0.63 0.23 0.14
## Cumulative Proportion 0.63 0.86 1.00
##
## Mean item complexity = 1.1
## Test of the hypothesis that 3 components are sufficient.
##
## The root mean square of the residuals (RMSR) is 0.07
## with the empirical chi square 9489.04 with prob < 0
##
## Fit based upon off diagonal values = 0.97
trip_distance, fare_amount,
duration, total_amount, tolls_amount’,
‘tip_amount’. Các biến này đều liên quan đến các khía cạnh chi phí và
khoảng cách của chuyến đi, bao gồm chi phí cơ bản, chi phí phụ (như phí
cầu đường và tiền tip), và các yếu tố thời gian và khoảng cách của
chuyến đi. Vì vậy, ta có thể đặt tên cho nhân tố thứ nhất này là “Chi
phí và Khoảng cách Chuyến đi”.Để kiểm tra thang đo nhân tố được rút ra được tạo thành từ các biến có phù hợp hay không ta sử dụng kiểm định Cronbach Alpha. Chẳng hạn để kiểm định xem nhân tố thứ nhất đặt tên là “Chi phí và khoảng cách” được cấu thành từ 6 biến có phù hợp không, ta sử dụng lệnh sau
psych::alpha(eda[, c('trip_distance', 'fare_amount', 'duration', 'total_amount', 'tolls_amount', 'tip_amount')])
## Number of categories should be increased in order to count frequencies.
##
## Reliability analysis
## Call: psych::alpha(x = eda[, c("trip_distance", "fare_amount", "duration",
## "total_amount", "tolls_amount", "tip_amount")])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.86 0.93 0.95 0.67 12 0.00046 8.1 6.6 0.63
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.86 0.86 0.87
## Duhachek 0.86 0.86 0.86
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## trip_distance 0.84 0.90 0.94 0.64 8.9 0.00055 0.037 0.62
## fare_amount 0.77 0.89 0.91 0.62 8.2 0.00080 0.031 0.63
## duration 0.82 0.92 0.95 0.68 10.8 0.00043 0.037 0.63
## total_amount 0.78 0.89 0.90 0.61 7.8 0.00093 0.032 0.58
## tolls_amount 0.89 0.93 0.96 0.73 13.8 0.00047 0.039 0.79
## tip_amount 0.88 0.94 0.95 0.76 15.8 0.00044 0.029 0.79
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## trip_distance 22699 0.92 0.92 0.92 0.91 2.91 3.7
## fare_amount 22699 0.98 0.96 0.98 0.96 12.90 10.5
## duration 22699 0.90 0.84 0.80 0.81 14.44 11.9
## total_amount 22699 0.98 0.98 1.00 0.97 16.12 12.9
## tolls_amount 22699 0.64 0.74 0.67 0.62 0.31 1.4
## tip_amount 22699 0.61 0.69 0.64 0.56 1.84 2.8
Một thang đo được xem là hợp lý nếu giá trị Cronbach Alpha > 0.7. Trong bộ dữ liệu này, giá trị Cronbach Alpha là ~0,86 nên có thể nói rằng thang đo “Chi phí và khoảng cách” được cấu hình thành từ 7 biến (items) là hợp lý.Tương tự, thực hiện kiểm tra thang đo cho 2 nhân tố còn lại.
psych::alpha(eda[, c('mta_tax', 'improvement_surcharge')])
##
## Reliability analysis
## Call: psych::alpha(x = eda[, c("mta_tax", "improvement_surcharge")])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.6 0.77 0.63 0.63 3.4 0.0031 0.4 0.025 0.63
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.59 0.6 0.61
## Duhachek 0.60 0.6 0.61
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## mta_tax 1.58 0.63 0.4 0.63 1.7 NA 0
## improvement_surcharge 0.25 0.63 0.4 0.63 1.7 NA 0
## med.r
## mta_tax 0.63
## improvement_surcharge 0.63
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## mta_tax 22699 0.97 0.9 0.72 0.63 0.5 0.039
## improvement_surcharge 22699 0.80 0.9 0.72 0.63 0.3 0.016
##
## Non missing response frequency for each item
## -0.5 -0.3 0 0.3 0.5 miss
## mta_tax 0 0 0 0 1 0
## improvement_surcharge 0 0 0 1 0 0
psych::alpha(eda[, c('passenger_count', 'extra', 'mta_tax', 'improvement_surcharge')])
## Number of categories should be increased in order to count frequencies.
## Warning in psych::alpha(eda[, c("passenger_count", "extra", "mta_tax", "improvement_surcharge")]): Some items were negatively correlated with the first principal component and probably
## should be reversed.
## To do this, run the function again with the 'check.keys=TRUE' option
## Some items ( passenger_count ) were negatively correlated with the first principal component and
## probably should be reversed.
## To do this, run the function again with the 'check.keys=TRUE' option
##
## Reliability analysis
## Call: psych::alpha(x = eda[, c("passenger_count", "extra", "mta_tax",
## "improvement_surcharge")])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## -0.0036 0.35 0.41 0.12 0.53 0.0057 0.69 0.34 0.02
##
## 95% confidence boundaries
## lower alpha upper
## Feldt -0.03 0 0.02
## Duhachek -0.01 0 0.01
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se
## passenger_count 2.2e-02 0.486 0.503 0.240 0.947 0.00211
## extra 7.4e-05 0.439 0.479 0.207 0.783 0.00079
## mta_tax -5.7e-03 0.032 0.023 0.011 0.033 0.00641
## improvement_surcharge -5.1e-03 0.037 0.027 0.013 0.039 0.00644
## var.r med.r
## passenger_count 0.11324 0.0508
## extra 0.13330 -0.0014
## mta_tax 0.00067 -0.0014
## improvement_surcharge 0.00109 -0.0063
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## passenger_count 22699 0.939 0.42 -0.0094 -0.0068 1.64 1.285
## extra 22699 0.335 0.47 0.0592 -0.0043 0.33 0.463
## mta_tax 22699 0.047 0.72 0.7196 0.0185 0.50 0.039
## improvement_surcharge 22699 0.042 0.72 0.7157 0.0306 0.30 0.016